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Abstract 

Stability of electro-hydrodynamic processes between ion-exchange membranes is inves- 
tigated. Solutions of the equilibrium problem are commonly described in the one- 
dimensional (ID) steady-state approximation. In the present work, a novel class of 
ID unsteady self-similar equilibrium solutions is developed, whose existence is supported 
by recent experiments. ID unsteady two-parametric family of self-similar equilibrium 
solutions and their stability are studied both asymptotically in small dimensionless Debye 
length, and numerically. The self-similar solutions and marginal stability curves obtained 
in both approaches are in a fair agreement with each other at intermediately large times 
for which the dimensionless distance between membranes is large. 
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I. INTRODUCTION 



Problems of electro-kinetics attract a great attention due to a rapid development of 
micro-, nano- and biotechnologies. Among numerous modern applications of electro- 
kinetics, note micro-pumps and biological cells, electro-polishing of mono- and poly- 
crystalline aluminum, growth of aluminum oxide layers for creating micro- and nano-scale 
regular structures, such as quantum dots and wires. 

Study of the space charge in the electric double-ion layer (EDL) is a fundamental 
problem of modern physics [I]-|S], first addressed by Helmholtz. A detailed description of 
the current state of the art is presented by Zaltzman and Rubinstein [9]. As distinct from 
the commonly accepted stability studies of the steady-state equilibria, the present work 
describes a novel class of unsteady transient equilibria predicted numerically (Demekhin 
et al. [To]) and observed experimentally (Yossifon and Chang ^\). Since the present 
approach is strongly based on the principal results developed for steady-state equilibria, 
they are discussed below along with basic experimental data. 

A typical experimental voltage-current (V-C) curve in Figure [1] gives a first insight 
into the system. The V-C curve depicts three characteristic regimes: the regime A of 
under-limiting currents proportional to a low potential drop; the regime B of limiting 
currents (the plateau region), and the regime C of over-limiting currents proportional to 
the potential drop as in the regime A. 

Theoretical description of ID steady-state equilibria. Early studies assume (see, e.g., 
|12] ) that outside of a small vicinity of the membrane surface y = 0, right away after EDL, 
there is an equilibrium diffusion layer which satisfies the condition of electro-neutrality. 
For ID steady-state equilibria homogeneous along membranes, this yields the following 
solution: 

^' -(y - 1) + 1, $ = In ■^'(^-^)+,^'- , AV = In jum = 2. (1) 



Jlim Jlim J Jlim J 

Here are ion concentrations, ^ y ^ 1 is the coordinate normal to the membrane 
surfaces, j, jum and AV^ have the same meaning as in Fig. [H The solution ([1]) is made 
dimensionless using the diffusion layer thickness, the bulk ion concentration and the 




Figure 1: Schematic V-C curve for ion-exchange membranes. A, B and C represent 
under-hmiting, hmiting and over-hmiting regimes, respectively; * depicts the threshold 
of the over-limiting regimes; C/ and Cn are the regions of regular and irregular current 
oscillations; j and jum are the current density at the bottom membrane and its limiting 
value; AV^ is the potential drop. 

thermodynamic potential as the characteristic scales. As follows from ([T]), the V-C curve 
obeys a linear Ohmic relationship for small j (region A in Fig. [1]) . As it was first pointed 
out by Levich [12], the solution ([1]) has no physical meaning for j ^ jum- 

Several works (Smyrl and Newman [Ij, Rubinstein and Shtilman [3], Grafov and Cher- 
nenko [13]) are focused to remedy this inconsistence. Rubinstein and Shtilman |i3j came 
up with ideas of the non-equilibrium nature of EDL and of the extended space charge 
(ESC) region which is much thicker than EDL. Besides this, the fundamental novelty of 
this study is the emergence of the diffusion layer near a permselective solid/liquid in- 
terface for the limiting and overlimiting currents. This becomes possible due to a lucky 
choice of a suitable model problem. The quiescent steady-state ID problem has a small 
parameter e (dimensionless Debye length) at the highest derivative of the Poisson equa- 
tion for the electric potential. The nonlinear problem is reduced to the Painleve equation 
of the second kind for the electric field and solved numerically. The work [3] gives a key 
to the understanding of the limiting regimes and provides the success of their further 
studies [I1]-[2I]. 

An effective asymptotic approach in small e (the decomposition method) yields the 
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following simple analytic solution of the ID quiescent steady-state problem |21]: 

c+ = = 0, F = AF 

= = —{y -Vm), F = 0, ym^y<l, (3) 

Jlim 

where the V-C characteristics and the thickness of the ESC region ym are 

AF = ^.l^yf/^ ym = l-^^, Jum = 2. (4) 

3 y Jlim J 

Here F = eV, AF = eA^. The solution ©-(jl]) may be additionally simplified for the 
near-critical values of ym and AF. These values are small at small current deviations 
(j ~ jlim), which in turn may be ordered in e. The relations ([2]) and ([3]) determine, to 
leading order in e, the outer solution of the problem in the space-charge and electro- 
neutral regions, respectively. In general, matching conditions should relate the outer 
solution with the inner solutions in the thin internal boundary layer at y = ym and in 
the double-ion boundary layer near y = 0. Here the terms "outer" and "inner" solutions 
have the usual meaning of the method of matched asymptotic expansions [23]. Instead, 
patching conditions are applied to the outer solution aX y = y^ and ?/ = 0. This yields 
a fair approximation of the problem solution in the outer regions applying the patching 
conditions to the outer solutions and ignoring inputs of the inner solutions in the boundary 
layers. This is qualitatively illustrated in Fig. [21 where the electric field E, the charge 
density p = c+ — and the concentrations versus coordinate y are depicted for both 
the asymptotic and exact solutions at small e. 

The solution describes extreme non-equilibrium over-limiting regimes, while 

([1]) is valid for the under-limiting regimes. A unified asymptotic description of EDL valid 
from under-limiting to extreme non-equilibrium over-limiting regimes is developed by 
Zaltzman and Rubinstein [5]. 

Experimental evidence of the electro- convective instability. With further increase of 
the potential drop between the membranes, the over-limiting currents eventually arise, 
and steady-state equilibria lose their applicability (transition from the regime B to C in 
Fig. [1]). In general, four physical mechanisms can be responsible for that phenomenon. 
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Figure 2: Schematic for limiting processes in a steady-state equilibrium: a) electric 
field E = d(^/dy\ b) volume charge density p= — c~\ c) ion concentrations c = c+ and 
c = c~. The dashed line is the outer solution in the small vicinities of y = (the double 
ion boundary layer) and y = ym (the boundary layer between the volume space charge 
region and the diffusion layer). The solid and dashed lines depict exact and asymptotic 
solutions for small e. 

Additional charge carriers due to electrolyte splitting and exaltation effects are assumed 
to be responsible for over-limiting currents in Belova et al. [23], Pismenskaya et al. [25] . 
Rayleigh-Benard convection (see review by Cross and Hohenberg pS]) can be also due 
to the Ohmic heating, hence, it can provide another physical mechanism responsible for 
the over-limiting regimes. However, for a relatively small distance between membranes, 
the gravity-convection effect is absent, since the Rayleigh number Ra is less than the 
critical value, i^a^ = 1708 (see experimental data in [23], [25]). Finally, Rubinstein, 
Staude and Kedem [27] found experimentally that the transition to the overlimiting 
regime C (see Fig. [T]) is accompanied by current oscillations which are regular at a 
small super-criticality and irregular at a large super-criticality (regimes Cj and Cu 
in Fig. [1]). The above-mentioned facts indirectly indicate to a correlation between the 
instability and over-limiting currents. The correlation between the electro-convective 
instability and the over-limiting regimes is demonstrated in Maletzki et al. [28] and 
Rubinstein et al. [29], where the over-limiting regimes are eliminated when instability has 
been artificially suppressed. The first direct experimental proof of the electro-convective 
instability, which arises with increasing potential drop between the bulk ion-selective 
membranes, is reported by Rubinstein et al. [30], who manage to show the existence of 
small vortices near the membrane surface. Yossifon and Chang also observe an array 
of electro-convective vortices arising under the effect of a slow AC electric field, while 
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Kim et al. [3T] observe non-equilibrium electroosmotic vortices. 

Model of electroconvection as a stability lost of steady-state equilibria. As it was 
mentioned above, the over-limiting currents are accompanied by an instability of the 
membrane system. The following two modes of electroconvection in strong electrolytes 
are distinguished: i) bulk electroconvection due to the volume electric forces [32]-[39]; and 
ii) common electroosmosis, either of the classical "first" kind |1Q]-[15] or of the "second" 
kind [in]^|13]- In the first kind electroosmosis, a slip velocity results from the tangential 
electric field applied to the space charge of a quasi-equilibrium EDL [9], while in the 
second kind regimes, slip velocity results from the field applied to the ESC. 

The stability theory of ID quiescent steady-state EDL is developed by Rubinstein 
and Zaltzman [16], [17] for the extreme nonequilibrium conditions ©-(jlj). First, the re- 
gion ^ y ^ L is subdivided into the space charge, ^ y < ym, and electro-neutral, 
ym< y ^L, regions. The latter consists of the diffusion layer region, ym< y <1, and 
the bulk neutral region, 1 < y ^L. For the over-limiting regimes, the current inhomo- 
geneity along the membrane leads to a convective motion of fluid with the tangential slip 
velocity, Um, which is determined from the solution in the ESC region, < y ^ ym- Then 
the concentration, pressure and velocity are determined in the electro-neutral region, 
ym< y ^L, using the slip velocity found above as the boundary value at y = ym- The 
ID steady-state equilibrium solution is found to be stable for under-limiting regimes and 
unstable for over-limiting regimes. The corresponding marginal stability curve (in the 
plane of the potential drop, A$, and wave number, a) has no upper branch, and the 
growth rate increases infinitely with the wave number. A regularized problem that takes 
into account the higher order terms in e and predicts the upper branch, is offered by Ru- 
binstein and Zaltzman [48j. Linear stability results for various asymptotic formulations 
along with numerics are presented by Rubinstein, Zaltzman and Lerman |49j . Instability 
theory of ID quiescent steady-state solution uniformly valid for both equilibrium and 
non-equilibrium conditions (not necessarily of the extreme type) is developed by Zaltz- 
man and Rubinstein [H]. They obtain that non-equilibrium electroosmotic instability 
occurs at the transition from the quasi-equilibrium to the non-equilibrium EDL regime. 
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Figure 3: Marginal stability curves at different e, numerics. The under-limiting regimes 
are stable (adopted from Demekhin et al. [TU]). 

Note the numerical simulations of the nonlinear unsteady problem which demonstrate, 
in qualitative agreement with experimental data, the formation of the pair vortices and 
their development to chaotic motion in the electro-neutral region ( [IHl |50] ) . 

New type of equilibria and their stability. Although the cells observed in Rayleigh- 
Benard convection and Benard-Marangoni t her mo-convection look like the cells in the 
electro-convective motion, the latter instability is much more complicated from both 
physical and mathematical points of view. In particular, electro-hydrodynamics of elec- 
trolytes between membranes manifests new types of equilibria and their instabilities de- 
veloped in micro- and nano-scales. 

The present numerical calculations evaluate the settling time as about several minutes, 
while experiments demonstrate that the characteristic time for the instability to manifest 
itself is about several seconds (see [TT], [21], [29] and the Subsection IIC). Thus, although 
the problem solution eventually reaches a steady-state equilibrium which becomes un- 
stable, an unsteady solution can lose its stability long before the relaxation to the steady 
state. These two unperturbed solutions, steady-state and self-similar unsteady ones, are 
rather correspond to different, gradual and instantaneous loading of the potential drop 
to the membrane system, respectively. 

The existence of this novel type of unsteady solutions has been recently predicted 
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numerically for small e at intermediately long times [10] and observed experimentally 
[TT] . In the present work, unsteady transient equilibria and their stability are investigated 
both asymptotically and numerically. The present asymptotic analysis is restricted by 
the limiting regimes and their transition to the over-limiting regimes. The results of the 
present asymptotic analysis are consistent with the exact numerical solution. 

The paper is structured as follows. In the next Section, the governing relations are 
presented and typical characteristic values of the system parameters are evaluated. The 
ID unsteady and ID self-similar problems are developed for intermediately long times. 
Numerical justification of closeness of the exact solution of the initial-boundary prob- 
lem to the solution of the self-similar problem is justified numerically for intermediately 
long times. In Section III, ID explicit solutions of the self-similar problem are obtained 
for small Debye length and compared with the numeric solutions for intermediately long 
times. A two-parametric family of self-similar unsteady equilibrium solutions is found 
and their stability is investigated. The justification of self-similarity is done both numeri- 
cally and asymptotically. In Section IV, the unsteady 2D problem is made dimensionless, 
and simplified equations are derived for small Debye length. Section V describes the 
linear stability of the ID self-similar solutions slowly varying with time. The resulting 
marginal stability curves are compared with those obtained by direct numerical simu- 
lations. Some of cumbersome calculations for electrostatic solution in the space charge 
region and for slip velocity are in Appendixes A and B, respectively. 



A. Governing equations and typical values of parameters 

Electro-convection in a binary electrolyte between semi-selective ion-exchange mem- 
branes is described by the equations for ion transport, the Poisson equation for the elec- 
tric potential and the Stokes equation for a creeping flow (this model was first introduced 



II. 



THE PHYSICAL MODEL AND 
ID SELF-SIMILAR PROBLEM 



in [9]): 
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V^i. = j{5-- 5+) , (6) 



- VP + /iV^U + F - 5+) V$ = 0, V ■ U = 0, (7) 

This system of dimensional equations is complemented by the proper initial and boundary 
conditions which are as follows 

5^=P, -^jr^ + ^ = 0, $ = 0, U = for y = 0, (8) 
RT oy dy 

Fc~ (9$ dc" ~ ~ ~ 

5+=P, -^77^ + ^ = 0, <f = Ar, U = for y = L, (9) 
RT ay dy 

c+ = 0- = Coo for t = 0. (10) 

Here c"*", c~ are the molar concentrations of cations and anions; U = {f/, V , W} is the 
fluid velocity; {x, 5} are the coordinates, y is normal to the membrane surface; $ is 
the electrical potential; e is the permittivity of the medium; P is the pressure; Coo is the 
initial ion concentration; F is Faraday's constant; R is the universal gas constant; T is 
an absolute temperature; jl is the dynamic viscosity of the fluid; D is the cationic and 
anionic diffusivity; p is the interface concentration. 

It is convenient to use the characteristic electric current j at y = which is related 
with the known total potential drop between membranes AV^ 

~ DF^ , 9$ ~ ~dc+ 

J = ^5+j^ + DF—, y = 0. 11 
RT oy dy 

Let us evaluate the parameters of the system. The bulk concentration of the aque- 
ous electrolytes varies in the range Coo = 1 -i- 10'^ mol/m^] the potential drop is about 
AV^ = 5 V; the absolute temperature can be taken as T = 300" K; the diffusivity 
is about D = 2 ■ 10^^ m? / s\ the distance between the electrodes L is of the order of 
0.5 ^ 1.5 mm; the concentration value p on the membrane surface must be much higher 
than Coo, and it is usually taken within the range from 5coo to lOcoo ( [3], [S], [SI])- The 
dimensional Debye layer thickness Ad is varied in the range from 0.5 to 15 nm depending 
on the concentration Coo- 
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B. ID unsteady problem at intermediately long times 

Let us consider an ID unsteady equilibrium state: 

d d 

Then the system flS])- ffTT]) turns into 

52$ p 



(14) 



Fr~ (9$ dr~ 

~c'=P. -^^ + ^ = 0, $ = for y = 0, (15) 
RT oy dy 

£+=P, -^^ + ^ = 0, «l> = AV^ for y = L, (16) 
RT oy dy 

= c~ = Coo for t = 0, (17) 
~ DF^ , a$ ~ ~dc+ 

3 = ^c^j^ + DF—, y = 0. 18 
RT oy Oy 

According to f ll3p -( TT5]) . at AV = 0, a neutral steady-state solution with a uniform 
concentration of ions, c"*" = c~ = Coo, instantaneously loses its equilibrium under the effect 
of the potential drop, AV^, applied at t = and kept constant at t > 0. As a result, a 
concentration polarization of ions is formed in thin boundary layers near the membrane 
surfaces, which expand with time. After several milliseconds of evolution, the diffusion 
layer thickness, S, will be much larger than the Debye length, Xd = {iRT /cooy^"^ / F, 
S{i) ^ Xd- On the other hand, there is a wide interval of time i (from seconds to 
minutes), such that 6{t) ^ L: 

Xd < < L. (19) 

We assume that the diffusion-layer thickness S{t) in the ID unsteady solution varies in 
time as 

5{i) = 2\fbi. (20) 
Then f|T9|) can be presented in the form 

(21) 
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Let us study the problem for intermediately long times which belong to the interval 
(^3). The spatial structure of the unsteady limiting regimes at any fixed time is qualita- 
tively similar to that of the steady-state equilibria depicted in Fig. |2j The space between 
the membranes can be divided into the following regions of inhomogeneity: (I) a thin 
EDL near the membrane surface y = 0; ESC region (II) separated from the diffusion layer 
(III) by a thin internal boundary layer (IV); a bulk neutral region (V). 

For unsteady solutions, the influence of the upper membrane y = Lis neglegible within 
the interval (12T]) . and the boundary conditions (flGl) should be revised. Subtracting the 
second equation (|T3|) for from the first equation (|T3l) for c"*", integrating the result from 
to some y within the bulk neutral region, and neglecting the higher order terms with 
respect to small S/L, the known linear voltage distribution in the bulk neutral region is 
obtained 

~ _ Irt ^ ^ ~ 

$(y) ~ const + ' y, y:^S. (22) 

Zr UCqq 

Here the constant and the electric current j are functions of time. The relation fl22|) 
predicts unbounded linear behavior of the potential at a large distance from the bottom 
membrane. This infers the following effective boundary condition at large y 

lim I !> - y\ = A$, (23) 

which, along with the following condition for the concentrations 

lim = lim c~ = Coo (24) 

y/S—^oo y/5—^oo 

substitute the EC's (|T6|) for the problem with the distance between the upper and bottom 
membranes that satisfies f[T^ . Note that Ohmic current is taken into account due to a 
cumulative effect at a large distance from the membrane. 

Introducing a similar effective voltage instead of the total voltage AV is useful also 
for the finite membrane problem both for the efficiency of calculations and comparison 
with the case of the distant upper membrane: 

A<l = Ay-4£^. (25) 
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In fact, the effective voltage A$ in the EDL and ESC regions does not depend on the 
distance between the membranes. 

The formulated problem does not contain a characteristic length, since neither the 
distance between the electrodes nor the double-ion length can give such a scale within the 
interval f l2T]) . This infers the solution self-similarity [52] with the diffusion layer thickness 
as the dynamic characteristic size ( justification of self-similarity in small e see in Section 
II. C). Let us transform the problem into a dimensionless form using the characteristic 
values of potential, ^o = RT/F, concentration, Coo, and time, to, by introducing the 
following new independent variables: 



This gives 





T = t, 


V -- 


y 






m ' 




dc+ 
dr 


dr] 


-( 

dr] \ 




, d^c+ 
' dr]"^ 


dc~ 
dr 


or] 


-( 

dr] \ 


' _d<i>\ 

V dr]j 


1 + 

' dr]'^ 








~-c\ 








= c 







(26) 



(27) 



where e = Xo/Sii). 

Assume now that for intermediately long times the solution, indeed, becomes self- 
similar (see Section II. C), i. e. 

„ rfc"*" d / , ^ dc~^ \ 
2r]— + — c+E + — = 0, 
dr] dr] \ dr] J 

dc' d f ^ dc-\ 
2r]— + — -c-E + — = 0, (28) 



dr] dr] \ dr] 



^dE 
dr] 



with the EC's 



dc~ 

c+ = p, c-E - — = 0, $ = for 77 = 0, (29) 
dr] 

c- ^1, c+ 1, $ - ^ A$ at 7/ ^ 00. (30) 
d^ ~ ~ 

Here E= ——, A$ = A$/$05 the dimensionless electric current at 77 = is determined by 
dr] 

j = c+E+^, r] = 0. (31) 
dr] 
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Figure 4: Self-similar charge p = — c (a) and concentration c"*" distribution (b) 
for different e. Curves correspond to the following values of electric potential drop: 1. 
A$ = 50; 2. A$ = 100; 3. A$ = 150; 4. A$ = 200. 

Solutions of the self-similar boundary problem fl28|) - fl3T|) are presented in Fig. H] for 
the charge distribution p{ri) = C^irj) — c~{ri) and positive ions concentration c^{ri) for dif- 
ferent values of the parameters. Voltage-current characteristics J versus A$ for different 
e are shown in Fig. |5] for both the exact numeric and asymptotic self-similar solutions. 

C. Numerical justification of the self— similar problem 

To justify self-similarity of the solution of (128|1 -( 13T]) for intermediately long times, we 
compare it with the exact numerical solution of the initially-boundary problem f[T^ - f[T51) . 
In this numeric simulations, the time t and the coordinate y are referred to / D and L. 
The concentrations and the potential $ are referred to Coo and $0; respectively, as it 
was done for the self-similar solution. The system f ll3p - (fT5]) turns into 



dt ^ dy^ dy) ~^ dy'^ ' ^'^'^^ 
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Figure 5: VC curves. Comparison of the self-similar asymptotics (dashed lines) and the 
exact numerics (solid lines) for several e. 




Figure 6: Distribution of the charge density p in space for several moments of time, 
AV^ = 100 and v = \d/L = 0.001. Solid line corresponds to the numerical solution of 
the unsteady problem Two-parametric self-similar solution of ( 12S]) -( PT]) with 

e and A$ parametrically changing with time is depicted by dashed lines. 



14 



(9$ dc~ 

=p, -c^— + — = 0, $ = for y = 0, (34) 
oy oy 

c+ =p, -c-— + — = 0, <I> = AV for y = l, (35) 
oy oy 

= c- = 1, for t = 0, (36) 

where v = \dIL, AV-Aiz/io, j^jL/{DFc^). 

The system fl32|) - fl371) is integrated numerically using Gear's method in time along 
with Galerkin's r-method digitization with respect to the space variable. The problem is 
described by two parameters, the voltage between membranes, AV, and analog of e in a 
finite size system, i/ ^ 1. The calculations are fulfilled for the limiting regimes, u and AV 
are varied in the ranges u = 0.00001-=-0.01 and AV = 50^200. It is found numerically that 
the solution reaches the steady-state equilibrium at ts= 0.1 ^0.2. Taking characteristic 
values D = 2 ■ 10^^ m^/s and L = 0.5 ^ 1.5 mm, this time is evaluated in a dimensional 
form as tg = 0.25 -i-4 min, supporting the assumption in Subsection AI. Below, we restrict 
ourselves by the regimes satisfying the inequality (!2T|) which can be generalized now as 

— <^t<^ts. (38) 



The two-parametric family of self-similar solutions of the problem f l28|) - fl3T|) is calcu- 
lated for a wide range of the parameters A$ and e and accumulated for further comparison 
with the exact numerical solution for the finite-length membrane geometry. For the com- 
parison of the self-similar solution with the exact numerical solution for the finite-length 
membrane geometry, the parameters A$ and e varying with time are calculated using 
relations 

.(t) = ^, J{t) = 2j{t)V~t, A$(t) = Ar-^. (39) 

Here the Ohmic voltage, j{t)/2, is excluded from the total voltage AV. These values 
are substituted into the accumulated family of self-similar solutions fl28l) - fl3Tl) obtained 
previously for a wide range of the two values A$ and e. 
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The charge density p = c"*" — c~ is presented in Fig. [6] for both exact and approxi- 
mate self-similar solutions. The unsteady solution is "adiabatically" sliding along the 
two-parametric self-similar solution, A$(t) and e{t). One can see a rather good corre- 
spondence between these two solutions for intermediately long times. 

This expands the applicability of the approach up to quite realistic values of the 
potential drop. Thus, the present asymptotic model is applied to the dimensionless 
potential drop < A$ < 200 (which corresponds to dimensional values -v- 5V). The 
asymptotic model describes the limiting regimes with the error value depending on e. For 
e = 0.01 and A$ = 40 the error is about 20%, while for A$ = 200 and the same e = 0.01 
it is about 3%. For e = 0.001 the error is much smaller, for example for A$ = 20 it is 
about 1% and for A$ = 200 less than 0.2%. 

For the self-similar solution the diffusion layer thickness should be inversely propor- 
tional to the electric current, 6(t) ~ l/j(i) ~ Vt- Typical evolution of l/j(t) presented 
in Fig. [7] exhibits three characteristic time intervals. In the interval //, for the interme- 
diately long times fl55]) . l/j(t) is proportional to \/i and the assumption fl2U]) is satisfied, 
while for short times, in the time interval J, as well as for long times, in the interval ///, 
the self-similarity is violated. 

III. ID SELF-SIMILAR SOLUTIONS FOR 

SMALL DIMENSIONLESS DEBYE LENGTH 

A. ID explicit solutions 

To find ID self-similar equilibrium solutions explicitly, we apply to (l25]) - fl3ip the asymp- 
totic decomposition method widely used for ID steady-state equilibria at small e, |14]- 
|21] . The reduced problem allows a relatively simple solution for the electric field, which 
consists of two branches on the left and right of the boundary of the space-charge region 
ri = rim- The value rim is determined by the solution of the problem. The resulting outer 
solution is weakly singular at = r^^, where the electric field gradient has a discontinuity. 
This procedure ignores the EDL near the membrane surface t] = (the first boundary 
condition c"^ = p in (129|) is rejected), but predicts the existence and location of the inter- 
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Figure 7: Typical evolution of l/j{t). I is the short-time region of influence of the initial 
data; II is the intermediately-long-time region of self-similarity, and III is the region of 
influence of the upper membrane y = l, {AV = 100 and z/ = A^/L = 0.001). 

nal boundary layer in the vicinity of i] = rim- The influence of these boundary layers on 
the limiting regimes is negligibly small for vanishing e. In particular, the potential drop 
within these boundary layers is much less by the order in e than that within the ESC re- 
gion. As in the steady-state case, the matching conditions of the outer solution with the 
inner ones within the boundary layers are substituted by the patching conditions. Such 
a description, in spite of its simplicity, well describes solutions of both equilibrium and 
stability problems as demonstrated below by their comparison with numerics. Within the 
adopted asymptotic approach, the interface elevation = 77^ is a new natural variable in- 
herent to the problem. It first arises explicitly in the equilibrium solution as the boundary 
between the extended space charge region and the region of the diffusion layer, in which 
the equihbrium states are described separately by the outer solution. Then the stability 
study employs perturbation of the interface elevation rjm along with the perturbations of 
all conventional variables such as concentrations, velocity etc. Hence, the system devel- 
opment is accompanied by the interface instability, which reflects the physical nature of 
the system. 

Following the decomposition method to leading order in e, the new scaled variables 
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F and F are introduced for the equilibrium problem: 



AF = £A$, H 



dF 
dr] 



(40) 



Substituting ( 140|) into ([28l)- (l3Ti) yields after some algebra: 

d'^r / 1 



dr]"^ 

d_ 
dfj 



r + -H^j H = 0, F = for r] = 0, 

r + ^//^ ^2, F ^ AF at 7] ^ oo, 
and the auxiliary condition for the electric current J related to F 



J 



(i?7 ' 



r] = 0. 



(41) 

(42) 
(43) 
(44) 

(45) 



To completely determine the outer solution in variables F, H and F, first the patching 
condition of the outer solutions is applied to the function (F + H^)H at 1] = !]^, ignoring 
the input of the internal boundary layer at ri = rim- Note also that the method of decom- 
position ignores the input of the boundary layer in r/ = into the outer region, and the 
boundary condition (29) for c"^ at r/ = is rejected. Then equations (H2|) and the first of 
equations fH3|) yield 

+ H = Q. (46) 

Equation ( HBj) yields the outer solution separately in the space-charge, ^ ^ 77^, and 
electro-neutral, rj ^ r^^, regions: 

1 



F + -H^ = for ^ < r/„ 



(47) 

if = for > r]^. (48) 
Then, using equations fH7|) - fH5]) . conditions of continuity for F, if, F and dV /dr] and 



the rest of equations fl4T|) - fl45|) lead to the following relations in the space-charge, and 
electro-neutral regions: 

3/2" 



r = J(r7-77j, F = AF 



1 



1 



Vm 



for ^ T] < T]„ 



(49) 
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r = 2 - 2 



eric^rj) 
erfc(?7m) 



F = AF = const for r] > rjm- 



(50) 



J 



Returning to the old variables, c"*", c , if and p, we obtain 
c+ = c" = 0, H = a/2 J (?7„ - r/), p = e 

erfc(?7) 



1 - 



erfc(r7ri 



2(^m - ^?) 
H = 0, p = for Tj > Tj. 



for ^ r] < rj„ 



Here, to the leading order in e, the volume charge in space-charge region has a singularity 
at ri = rim- To overcome this inhomogeneity of the outer solution f H^ - flSUj) in the vicinity 
of ri = rim, the inner solution should be developed. 

According to this solution the potential drop in space-charge region is approximately 
equal to the total voltage between the membranes: 



AF = ^ eAV. 



(51) 



The values AF and J are interrelated: 



AF 




^ 'Im 1 



(52) 



with the relation between r/^ and current J 

4 exp (— r]. 



J 



(53) 



a/tt erfc(r7m) 

Relations ( 152|) and (153!) determine implicitly the V-C curve and the boundary of the ESP 
region rjm versus potential drop AF 



Jerfc 



/9AF2 
V 8J 



2\ 1/3' 



Jlirn exp 



/9Air2\2/3- 

I 8J 



^^AF^ - ^mexp(-r7; 



^iim = (54) 
^/^T 



(55) 



32 erfc {r]m) 

In the subsequent stability study of the self-similar regimes, the scaled value of the 

potential drop is assumed to be small: 



AF = £A$ ^ sAV < 1. 
In that limit relations flMl)-fl55D are reduced as follows: 

2n. 



(56) 



J ~ Jlim I 1 ~l~ 



TT 



Vm ^ ^ {J - Jl 



lim ) 



(57) 
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Figure 8: VC curves. Solid line relates to (EH), while dashed line - to the simplified 
version flSSjl . 

{AF)' = ^{J-Juj'j. (58) 

Evidently, f l57p - (|H5]) correspond to near-critical currents. On the other hand the analysis 
of the inner solution at t] = rjm, and of the next order in e terms in the outer solution 
at T] > rjm yields that the approximate solution (H9l) - (l50i) homogeneously satisfies the 
condition flHTl) of the potential drop domination in space-charge region if 

£ln£-i<AF<l. (59) 
Then relations fl56l) - fl59p gives for a small potential drop: 

e\ne-'<^ri'J^r^iJ~JumY^^r^AF<^l, Jum = 0{{AFf). (60) 
For a small potential drop, relations ( H^ and f l5U]l yield within space-charge region: 

c+ = c" = 0, F ^AF F ~ (AF)2/3, H ~ (AF)l/^ p ~ (AF)-^/^, (61) 
and electro-neutral region: 

c+ ~ ~ (AF)°, F AF, F ~ (AF)°, H = 0, p = 0. (62) 
while according to equation (IBUj) . the boundary between these two regions is given by 

(e In £-^)2/3 < ~ ( AF)^/^ ^ ^_ ^gg) 
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Figure 9: Negative ion concentration (a), c , and electric field (b), if, vs y for several 
moments of time. Comparison of numerics f lT5]) -f lTH]) (solid lines) with self-similar asymp- 
totics (in])-(ISnD (dashed lines) for (a) A\/ = 50, z/ = 0.001: 1. t = 0.0001; 2. t = 0.001; 3. 
t = 0.004; 4. t = 0.015; 5. t = 0.5; and (b) AV = 100, = 0.0005: 1. t = 0.001; 2. t = 0.003; 
3. t = 0.03. 

B. Comparison of ID explicit asymptotic and numeric solutions 

In Fig. O the asymptotic solution for V-C curves is presented as a dependence J 
versus A$ (dashed lines) and compared with those (solid lines) obtained numerically in 
Subsection II. C. The comparison demonstrates a fair correspondence between the self- 
similar asymptotics and exact numerical solution except the region of under-limiting 
currents. The difference tends to zero as £ — e.g. at e = 0.001 the exact solution and 
asymptotics coincide with graphic accuracy. This justifies the neglect of the EDL inputs 
into the limiting solutions accepted in the present study. 

In Fig. |8] the universal voltage-current characteristic described by ( 15^ is plotted by 
solid lines, while simplified version of the VC curve f lSS]) is plotted by dashed lines. It is 
seen that the difference between these two curves is extremely small for low values of AF 
and does not exceed 4% at AF = 1. 

In Fig. in] the distributions of the negative ion concentration c~ and electric field 
H = eE obtained from the numerical solution of ( I32|) - (l36|) are shown for several t. In 
turn, in Fig. |9]d self-similar distributions described by ( H9|) -( [50|) are given. Except for 
short and long times, a rather good correspondence is observed of the both approaches 
at intermediately long times. In accordance to fl63|) . the value of rjrn, indeed, is a small 
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2 

Figure 10: Comparison of the V-C characteristics for self-similar solution and numerics, 
AV = 50, z/ = 0.00005. Numerics (lT51)- (fT5]) is shown by squares and the universal self- 
similar VC curve ( 154|) - by solid line. The region of the solution self-similarity ti< t <t2, 
ti ^ i^^/4, t2 '^ts = L'^/{4:D). Arrows depict the time growth along the voltage-current 
curve. 

value for parameters adopted in Fig. |9l where rjm is determined by the condition that 
r] > rjm correspond to positive values of c~. 

The self-similar character of the solution for intermediately long times is illustrated 
in Fig. [TOl In this calculation = 50, z/ = 0.00005 are kept constant during numerical 
integration of the system flH21) - fl5Bl) . Also, the potential drop between the membrane is 
kept constant, its fraction A$ = AV — j{t)/2 is changing in time as well as J. It is clearly 
seen that the numerical solution follows the self-similar equilibrium. This fitting persists 
up to the moment t = t2, when the influence of the upper membrane becomes significant. 

The data for several AV are gathered in Fig. [TTk ) and [TTb) for u = 0.001 and 
1^ = 0.0005. The data for short times t<ti and for long times t>t2 are discarded from 
the consideration and only data for intermediately long times are kept. The numerical 
points for all AV and u shrink into the universal VC curve, J{AF), described by (1541) . 
It is seen that the fitting is better for smaller i> and AV. 
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Figure 11: Numerical points for several values of the potential drop between the mem- 
branes, Al^, shrink into the universal VC curve J(AF) given by f l54p : a) z/ = 0.001, b) 
zy = 0.0005. 

IV. 2D UNSTEADY PROBLEM FOR SMALL DI- 
MENSIONLESS DEBYE LENGTH 

A. Governing relations. Dimensionless variables 

Aiming at a stability study of the ID equilibrium self-similar solutions, the 2D unsteady 
problem is rewritten using self-similar variables. First, the characteristic length scales, 
/o and 5, tangential and normal to the membrane surface are introduced together with 
the characteristic scales for stream function, velocity components, pressure and electric 
potential as follows: 



Pq = $0 = — ^. (64) 



The characteristic time scale, to, is also introduced as an arbitrary moderately long time 
at which the equilibrium solution is already settled on the self-similar regime and yet 
does not leave it; 5{t) is the diffusion layer thickness, while a is the characteristic wave 
number. Preserving the same notations for dimensionless stream function, velocities and 
pressure as for their dimensional values, the following independent self-similar variables 
are introduced for the 2D problem: 



t 



X 



y 



(65) 
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Then the time-dependent wave number and Debye length are introduced for further 
convenience 

«(r)=5(t>, e{T) = ^ (66) 

d{t) 

Functions q;(t), £:(t), characterizing the time variation of the equihbrium solution, are 
assumed to be slowly varied in the fast-time scale of perturbations. 
Nonlinear 2D problem may be written in the following form: 

DK 



Dt 



V ■ (pV$) + V'K, (67) 



^ = V- (irV*) + VV, (68) 
JJt 

e^v^cj) = -p^ (69) 

V^^ = -^V X (pV'$), (70) 



Here 



V - (a— — ^ V- (— — - 4r— - 2ri— + (u— + V— 

K = + c~ , p = c'^ — c~ as previously in ID problem, x = e^l/{p,D) describes physical 
properties of the liquid and electrolyte; is the stream function: 

^ dr] ' ^ dx ^'^ ^ 

Assuming the problem self-similarity, the following boundary conditions are accepted: 

9$ dK dp d<iJ 
$ = 0, iK-p)-- — + ^=0, K + p = 2p, ^ = ^=0 for, = 0, (72) 

Jri 

$ - Y ^ P ^ 0, K ^2, ^ -> 0, ^ ^ at T]^ oo, (73) 



where J is the electric current f l3T|l . 

B. Decomposition method for 2D problem. Ad— hoc model 

Extending the decomposition method for the limiting and over-limiting regimes from ID 
case to 2D case, the following auxiliary variables scaled in e are introduced: 
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First the system fl67|l -( 170|l is rewritten in the variables substituting p = — eV^F: 

DK 

-V ■ [(V'F)VF] + V^K, (75) 



Dt 

^,D{V^F) ^ ^^^^^ ^ ^2^4^^ ^^g^ 



Dt 



V^^ = — V X [{V^F)VF], (77) 



and then substituting K = T + ^H"^ 



= -V ■ [(V^F) VF] + + ^H' ) , (7J 

1 



e'^^^^p- = -V ■ (^r + ^//^J VF I + e'V^F. (79) 

The following ad-hoc model is adopted for the space-charge and electro-neutral re- 
gions, respectively: 

K = T + ^H^ = for 0<:r]<r]m, (80) 

F = AF = const {H = 0, p = 0) for r]^ > r]. (81) 

Substituting f lHOjl and f lSTl) into fl75l) - fl77|) and fl78|) - fl79|) . respectively, for the space-charge 
and electro-neutral regions yield taking into account the boundary conditions (172]) - (1751) : 



T+Ih^ = 0, V-[(V^F)VF] = 0, V^^ = ^V'x[(V2F)V'F], ioT 0^r]<7]^, (82) 
2 

DT 



F = AF = const, —- = V^T, V^^ = for r]^ > rj. (83) 
Dt 

The systems (l82l) and (l83l) should satisfy the boundary conditions adopted from fl72|) - (173i) : 

F = 0, ^ = — = for r/ = 0, (84) 
or] 

r ^ 2, ^ ^ 0, ^ at ?7 ^ oo. (85) 

OT] 

As in ID model, the method of decomposition ignores the input of the boundary layer 
in the near-bottom region into the outer region, and rejects the boundary conditions for 
c+ = (ir + p)/2 at ?7 = in ([72D. Additionally, the problems §2l> and (jHS]) should be 
complemented by the proper patching conditions which provide continuity of the solutions 
of the systems (182|) and (l83ll at the patching point t] = rj^: 

r_ = r+ = o, F_ = F+ = AF, v„F_ = o, v„r_ = Vnr+, v^^_ = v^^+, (86) 

25 



where V„ is the gradient projection to the direction normal to the interface rj = rjm'-, 
AF = const; k = 0, 1, 2, 3 due to the fourth order in rj of the differential equations for 
\E'; subscripts — and + denote the values of the corresponding variables on the left and 
right of 1] = rjm, respectively (for details see [22]). 

The resulting non-linear model consists of a small parameter e, but the linearized 
problem is further simplified by an additional rescaling in the electroneutral region that 
excludes e from the linear problem. In the next section the stability problem obtained 
by linearization about the self-similar ID solution is considered. As is clear from the 
above analysis of the ID self-similar problem, the resulting unsteady problem depends on 
two small dimensionless parameters: Debye length and total potential drop, and further 
modeling is restricted by small values of the potential drop. Finally, the results for the 
asymptotic model are compared with direct numerical simulations for the exact problem. 

V. STABILITY OF SELF-SIMILAR SOLUTIONS 

The developed ID solution can lose stability, and a new electro-convective regime bifur- 
cates at some values of parameters. Let us restrict ourselves with the marginal stability 
problem, 

d 

or 

and impose infinitesimal perturbations on the self-similar ID solution. Using the fact 
that the coefficients of equations fl82l) - fl83|) are slowly changing at long times: 



R = R(n) + R{ri) exp {ix) , (87) 

where R stands for any physical variables, e.g. for F, H, F, rjm, J or \E'. Bars and 
hats denote the unperturbed and perturbed variables, respectively (below bars for the 
unperturbed variables are dropped wherever possible with no confusion). 

A. Space charge region, ^ r] < rj^. 

Electrostatic problem. Substituting the relations and fl57|) into flS^ and linearizing 
with respect to small perturbations f , H, F, fjm, J, yields the following system in the 
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space charge region, ^ t] < rjm- 

f + = 0, (88) 

dr] 

— ^a"^ (h— + 2—f] - (89) 
dr]'^ \ drj drj J 

Linearizing the BC's (Ell), dSS]) and shifting them to the undisturbed boundary t] = 1]^, 

yields the boundary conditions for F and F: 

^ = 0, ^ = J for r/ = 0, 

dr] (90) 

F = 0, f + Jfjm = for T] = r]m. 
Four conditions (190|1 serve the boundary conditions for the third order system (I88l) - (l89|) 
and determine the disturbed electric current J. After eliminating F, our system (IHSjl - fl^U]) 
turns into the boundary problem for F: 

d^ ( jjdF\ ^2 ( ^dF ^dHp\_^ 
d-rf \ dr] J \ dr] dr] j 

F = 0, ^{ H^] =-J for r] = 0, 

dv\ dr]) ^g2) 
dp 

F = 0, H— = Jf]m for r] = r]m. 
dr] 

Our knowledge of the solution is complemented by the following expression for the charge 
density perturbation, p: 

Hydrodynamic problem. The linearized hydrodynamic problem in the region ^ r/ < 
?7m, subjected to the boundary conditions (15^ at r] = r]^ and at r] = 0, is as follows: 

^ _ 2^2^ + c,4^ - - ^F - a^Hp] (94) 

dr]'^ dr]"^ I drj^ drj^ j 

d^ 

^ = — = for 77 = 0, (95) 
dr] 

d^^ d^^ , , ^ 

— - + 2a— + = 0, —r-'ic?—-2o?^ = {) for r] = r]^. (96) 
dr]'^ dr] dr]-^ dr] 

Solution of the problem ( 1M1) -( 196|) determines the perturbations of the shp velocity 

d^ 

Ura = Vm = for r] = r]m. (97) 

drj 
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Figure 12: Components of scaled slip velocities u and v. 

Solutions of the electrostatic and hydro dynamic problems in the space charge region 
are presented in Appendixes A and B. In particular, the scaled slip velocities, u and v 
are depicted in Fig. [12] 

= -laCu (r) , V^ = Cv (r) . (98) 

where C* is a constant factor (see Appendix B). For small r tangential component of 
velocity coincides with obtained in [47] . The values of the components of the scaled slip 
velocities at rj = rjm will be used as boundary conditions for the hydrodynamic problem 
in the electro-neutral region rj > r]m- 

B. Electro— neutral region, r] > rjm- 

Linearizing the problem (155]) . (^Si), (JHS]) about the ID self-similar solution yields the 
eigenvalue stability problem in the electro-neutral region rj > rjm- 
Electrostatic problem. 

^ + 2r7^ - a^r = -ia^ (r^) e"--" , (99) 
df]^ drj 

__2„^_ + „.*=0, (100) 
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The problem is completed by the following boundary conditions: 



A dV dT 

T + Jr]m = 0, -1- = -r 
arj drj 



+ 2Jr]mfim, '^=lVm, = V = Vm, (101) 

ri=rim-0 



r ^ 0, ^ ^ at r] ^oo, (102) 

where the slip velocity Um, Vm is taken from the solution of the problem in the region 
^ < at ?7 = ?7„. 

Hydrodynamic problem. Notice that the hydrodynamical part of the problem can be 
solved separately and independently from the rest part of the problem. Solution in the 
space-charge region is presented in Appendix B, Eqs. (11321) . (11351) — (I138P (for details see 
also [22]). The velocity distribution at r/m < V < oo can be found by solving (llOOp with 
the boundary conditions (110 II) and (11021) . Solution is found using (1^ and 

^ = iCifj, 

[v + iu- v)^] e« + ^piO , 0<v^Vm, (103) 
[v + {u- v)^] T] > r]m, 

where ^ = airim—rf), the values (7, m, f , ipp are calculating according to Eqs. (11321) . (I139p — 

dm]). 

Eventually, the stream-lines determined as the level curves of the imaginary part of 
the perturbed stream function can be built up 

\1> = iiiji) exp(2x), ^ = iC%lj, (104) 

Tie (*) = -2Ci){T]) sin(x) = c(r). (105) 

Note that the stream-function lags behind the charge density by 7r/2, while the charge 
density p is in a counter phase with the interface elevation r]m- Since the function \l/ 
has a single extremum point and vanishes at and infinity, the stream-lines have closed 
trajectories corresponding to 2D electro-convective rolls. The calculation results for 
A$ = 100, £ = 0.001 are presented in Fig. [T3]a)-f). For increasing a, the centers of the 
rolls arriving from the infinity are crossing the boundary of the space charge region, and 
then, at a — )■ oo, are approaching the boundary rj = rjm from below. 
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Figure 13: Stream-lines for A$ = 100, e = 0.001 and different values of a. 
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C. Results for the linear stability analysis. 

The final step of our analysis is to find T in the electro-neutral region, rj > rj^n, and solve 
the stability problem. Boundary conditions for F at = 77^ are 

r = J, — = - {3aa [r) + 2r]m) J for r] = rfm- 

a drj a 

Further it will be suitable to normalize F to GJ /a, 

^ + 2r/^ - «2p ^ ^ ^cr^^Wr. u^ciiu-v) (rim - v)) e""^"^' , (106) 
dr]'^ drj ae^ 

f = 1, — = -3aa - 2?7^ for r] = r]m, (107) 
drj 

f ^ at 77 00. (108) 

The present asymptotic analysis in small e for the limiting regimes has no restrictions 
for the value of a. The system fll06l) - fll07p is solved numerically by the shooting method, 
and the calculated marginal stability curves are presented in Fig. O In this picture 
results of the numerical analysis without any restriction on a and e are presented [TOj . 
The error of the present asymptotics can be significant near the nose of the marginal 
stability curves, but small far from it. This occurs since the nose of the marginal curve is 
located near the transition region from the under-limiting to limiting regime (see Fig. [3]), 
where our asymptotics is not accurate. 

Since the dimensionless Debye thickness and wave number are slow functions of time at 
moderately long times, we obtain, using for clearness the dimensional variables (denoted 
by tilde): 



a = 



Note that the coefficients of the eigenvalue problem fll06l) - fll07p are slowly varied with 
time. Note also that in spite of their slow changing, the parametric dependence of a and e 
on time significantly changes the interpretation of the marginal stability curves. Taking 
this into account and following Shtemler (1981) [53], we present our neutral curves in 
coordinates \je and a, taking A$ as a constant parameter. The instantaneous values 
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Figure 14: Marginal stability curves for self-similar asymptotic in e solution fll06p - fll07p 
(dashed lines) and exact numerical solution (solid lines) [TU] . 



oil/e and a depend on time implicitly, a ~ and 1 ~ with the constant 

ratio 

Hence, any straight line which starts in the origin of the plane (1/^, a) with inclination 
k 

a = k- (109) 

€ 

characterizes wave number a constant with respect to time. If time is increasing, the 
instantaneous values of l/e and a lying on this line are increasing as For large 

enough k, the whole straight line is in the stable region. Then, if we decrease the tangent 
k, at some k = ko the line fll09p will be tangential to the neutral curve (see Fig. [151 upper 
left corner). The dependences of ko, = a(A;o) and I/eq = l/e^ko) as functions of A$ 
are given in Table [TJ 

With further decrease of k < k^, the straight line crosses the neutral curve in two 
points, 1 and 1', and the region along the line can be divided into three regions: I - the 
region between the origin and point 1; II - the region within the neutral curve, between 
points 1 and 1'; III - the region lying above the upper branch of the neutral curve. The 
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«o 


l/eo 


Q(A$) 


i?(A$) 


15 


0.02 


4.38 


218.8 


0.14 


1.83 


20 


0.06 


3.36 


57.9 


0.24 


0.96 


30 


0.104 


2.69 


25.85 


0.33 


0.526 


45 


0.14 


2.43 


17.62 


0.4 


0.321 


70 


0.15 


2.43 


14.17 


0.43 


0.194 


100 


0.17 


2.43 


16.08 


0.45 


0.134 



Table 1: Critical parameters which separate stable and unstable regions (tangent ko, wave 
number ao and l/eo); and Q and R are functions of A$. 
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initial point {l/e,a) belonging to region I corresponds to decaying perturbations, but 
with increasing time the instantaneous values of {l/e,a) intersect the lower branch of the 
neutral curve in point 1 and enter the instability region. For {l/e,a) in region II, the 
perturbations are increasing until the instantaneous values [l/e^a) intersect the upper 
branch of the neutral curve; after this instant moment they begin to decay. For [l/e^a) 
in region III, the perturbations decrease monotonically at all times. 

It is possible to show that for k < the straight lines always leave the instability 
region. Indeed, from ( 152|) . (I53|) in combination with ( 157|) . ( 158|) the asymptotics for the 
V-C curve as e — t- can be obtained 

just from this relation 

r = «r7^ ~ (hl^^^ a^l\ (110) 

Now the upper branch of the marginal stability curves in short-wave approximation, 
a — )■ oo, can be readily obtained. In particular, we can claim that 

/ IX 2/3 

a~Q(A$) - , (111) 



where the function Q(A<I>) is calculated numerically and presented in Tabled] 
In order to evaluate the lower marginal stability branch, it is supposed that 

a ~ /?(A$). (112) 

At e — )■ the relation becomes exact. The dependence i?(A$) is tabulated in Table [TJ 

The relations fillip and f lll2p allow, in particular, to estimate the time gap At in 
which the perturbation with k < ko belongs to the instability region, between points 1 
and 1', in Fig. [151 

,~ ~, ~ k'Q'-R^ 
At = t' -t — . — . 

This relation becomes exact at k = aXn ^ 0. 
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VI. CONCLUSIONS AND DISCUSSION 



Electro-convective processes in an electrolyte solution between a semi-selective ion- 
exchange membranes have been investigated numerically and asymptotically in the limit 
of small Debye length e. First, a simplified system of equations has been derived from the 
full system of equations that describes the ion transport and contains the Poisson and 
Stokes equations. Asymptotic expansions in small dimensionless Debye length are ap- 
plied to both equilibrium and stability problems for the limiting regimes. Only the outer 
solution has been considered, ignoring the input of the surface and internal boundary lay- 
ers. A novel class of the ID unsteady self-similar solutions is found. A linear stability of 
these solutions and transition to electro-convection is considered. The marginal stability 
curves are constructed both numerically and asymptotically. 

The results of stability modelling in small e are in fair agreement with those obtained 
by direct numerical simulations of the entire problem everywhere except for the vicinity 
of the critical point of the marginal stability curves. The possible reason of that deviation 
is that the critical point of the marginal curves is located near the transition region from 
the under-limiting to limiting regimes, while our asymptotics is valid only for the limiting 
regimes. 

ACKNOWLEDGMENTS 

E.A.D. is grateful to I.Rubinstein for his stimulating questions and remarks and to V. 
V. Nikonenko for numerous discussions and significant help in understanding the physical 
pattern the process. 



35 



References 



[1] W. H. Smyrl and J. Newman, "Double layer structure at the limiting current," 
Trans. Faraday Soc. 63, 207 (1967). 

[2] R. P. Buck, "Steady state space charge effects in symmetric cells with concen- 
tration polarized electrodes," J. Electroanal. Chem. Interf. Electrochem. 46, 1 
(1973). 

[3] I. Rubinstein and L. Shtilman, "Voltage against current curves of cation exchange 
membranes," J. Chem. Soc. Faraday Trans. II 75, 231 (1979). 

[4] V. V. Nikonenko, V. I. Zabolotsky and N. P. Gnusin, "Electric transport of ions 
through diffusion layers with impaired electroneutrality," Sov. Elektrochem. 25, 
301 (1989). 

[5] A. V. Listovnichy, "Passage of currents higher than the limiting one through the 
electrode-electolyte solution system", Sov. Electrochem. 51, 1651 (1989). 

[6] J. A. Manzanares, W. D. Murphy, S. Mafe and H. Reiss, "Numerical simulation 
of the nonequilibrium diffuse double layer in ion-exchange membranes", J. Phys. 
Chem. 97, 8524 (1993). 

[7] V. A. Babeshko, V. I. Zabolotsky, N. M. Korzhenko, R. R. Seidov, and M. A.- 
Kh. Urtenov, "The theory of steady-state transport of binary electrolyte in a 
one-dimensional case", [in Russian] Dokl. Akad. Nauk 361, 41 (1998). 

[8] K. T. Chu and M. Z. Bazant, "Electrochemical thin films at and above the 
classical limiting current," SIAM J. Appl. Math. 65, 1485 (2005). 

[9] B. Zaltzman and I. Rubinstein, "Electroosmotic slip and electroconvective insta- 
bility", J. Fluid Mech. 579, 173 (2007). 

[10] E. A. Demekhin, E. M. Shapar, and V.V. Lapchenko, "Initiation of electrocon- 
vection in semipermeable electric membranes", Doklady Physics 53, 450 (2008). 



36 



[11] G. Yossifon and H.-C. Chang, "Selection of nonequilibrium overlimiting currents: 
universal depletion layer formation dynamics and vortex instability" , Pliys. Rev. 
Lett. 101, 254501 (2008). 

[12] L. G. Levich. Physicochemical Hydodynamics, Prentice-Hall, New York, 1962. 

[13] B. M. Grafov and A. A. Chernenko, "Theory of passage of a constant current 
through a solution of a binary electrolyte" , [in Russian] Dokl. Akad. Nauk SSSR, 
146, 135 (1962). 

[14] M. A.-Kh. Urtenov, E. V. Kirillova, N. M. Seidova, and V.V. Nikonenko, "De- 
coupling of the Nernst-Plank and Poisson equations. Applications to a membrane 
system at overlimiting currents", J. Phys. Chem. B 111, 14208 (2007). 

[15] V. A. Babeshko, V. I. Zabolotsky, E. V. Kirillova, and M. A.-Kh. Urtenov, 
"Decomposition of Nernst-Planck-Poisson equations," [in Russian] Dokl. Akad. 
Nauk 344, 485 (1995). 

[16] V. A. Babeshko, V. I. Zabolotsky, N. M. Korzhenko, R. R. Seidov, and M. A.- 
Kh. Urtenov, "Theory of stationary transport of a binary electrolite in one- 
dimensional case", [in Russian] Rus. J. Electrochem. 8, 863 (1997). 

[17] V. A. Babeshko, V. I. Zabolotsky, N. M. Korzhenko, R. R. Seidov, and M. A.- 
Kh. Urtenov, "Theory of stationary transport of a binary electrolite in one- 
dimensional case. Numerical analysis," [in Russian] Dokl. Akad. Nauk 355, 488 
(1997). 

[18] V. A. Babeshko, V. I. Zabolotsky, R. R. Seidov, and M. A.-Kh. Urtenov, 
" Decomposition equations for the stationary transport of electrolite in one- 
dimensional case", [in Russian], Rus. J. Electrochem. 33, 855 (1997). 

[19] M. A.-Kh. Urtenov and R. R. Seidov. Mathematical models for the electrical 
membranes for water purification, [in Russian] KubSU, Krasnodar, 2000. 



37 



[20] M. A.-Kh. Urtenov, "Mathematical models for the electrical membranes for water 
purification", [in Russian] Doctor Thesis, KubSU, Krasnodar, 2001. 

[21] M. A.-Kh. Urtenov. Boundary problems for the electrical membranes for water 
purification, [in Russian] KubSU, Krasnodar, 1998. 

[22] S. V. Polyanskikh, 'Electrohydrodynamic stability of some microti ows with con- 
centration polarization'", [in Russian] Ph. D. Thesis, Moscow State University, 
Moscow, 2010. 

[23] M. Van Dyke. Perturbation Methods in Fluid Mechanics, Academic Press, New 
York, 1964. 

[24] E.I. Belova, G.Yu. Lopatkova, N.D. Pismenskaya, V. V. Nikonenko, C. Larchet, 
and G. Pourcelly, "Effect of anion-exchange membrane surface properties on 
mechanisms of overlimiting mass transfer", J. Phys. Chem. B 110, 13458 (2006). 

[25] N.D. Pismenskaya, V. V. Nikonenko, E. I. Belova, G. Yu. Lopatkova, Ph. Sistat, 
G. Pourcelly, and K. Larshe, "Coupled convection of solution near the surface of 
ion-exchange membranes in intensive current regimes", Rus. J. Electrochem. 43, 
307 (2007). 

[26] M.C. Cross and P. G. Hohenberg, "Pattern fomation outside of equilibrium". 
Rev. Modern Physics 65, 3, 851 (1993). 

[27] I. Rubinstein, E. Staude, and O. Kedem, "Role of the membrane surface in 
concentration polarization at ion-exchange membrane". Desalination 69, 101 
(1988). 

[28] F. Maletzki, H.W. Rossler, and E. Staude, "Ion transfer across electrodialysis 
membranes in the overlimiting current range: stationary voltage-current charac- 
teristics and current nouse power spectra under different condition of free con- 
vection", J. Membr. Sci. 71, 105 (1992). 



38 



[29] I. Rubinstein, B. Zaltzman, J. Pretz, and C. Linder, "Experimental verification 
of tlie electroosmotic meclianism of overlimiting conductance through a cation 
exchange electrodialysis membrane", Rus. J. Electrochem. 38, 853 (2004). 

[30] S. M. Rubinstein, G. Manukyan, A. Staicu, I. Rubinsten, B. Zahzman, R. G. H. 
Lammertink, F. Mugele and M. Wessling, "Direct observation of nonequihbrium 
electroosmotic instability", Phys. Rev. Lett. 101, 236101 (2008). 

[31] S. J. Kim, Y.-C. Wang, J. H. Lee, H. Jang and J. Han, "Concentration polar- 
ization and nonlinear electrokinetic flow near a nanofiuidic channel" , Phys. Rev. 
Lett. 99, 044501.1 (2007). 

[32] A. P. Grigin, "The convective coulombic instability of binary electrolytes in cells 
with plane-parallel electrodes", Sov. Electrochem. 21, 48 (1985). 

[33] A. P. Grigin, "Coulomb convection in electrochemical systems", Sov. Elec- 
trochem. 28, 247 (1992). 

[34] I. Rubinstein, T. Zaltzman, and B. Zaltzman, "Electroconvection in a layer and 
in a loop", Phys. Fluids 7, 1467 (1995). 

[35] R. Bruinsma and S. Alexander, "Theory of electrohydrodynamic instabilities in 
electrolytic cells", J. Chem. Phys. 92, 3074 (1990). 

[36] R. S. Alexandrov, A. P. Grigin, and A. P. Davydov, "Numerical study of electro- 
convective instability of binary electrolyte in a cell with plane parallel electrodes" , 
Russian J. Electrochem. 38, 1097 (2004). 

[37] J. C. Baygents and F. Baldessari, "Electrohydrodynamic instability in a thin 
fluid layer with an electrical conductivity gradient", Phys. Fluids 10, 301 (1998). 

[38] M. E. Buchanan and D. A. Saville, "Electrohydrodynamic stability in electro- 
chemical systems", in Proceedings of APS 53rd Annual Meeting, Nov. 19-21, 
Washington, DC (2000). 



39 



[39] I. Lerman, 1. Rubinstein and B. Zaltzman, "Absence of bulk electroconvective 
instability in concentration polarization", Phys. Rev. E., 71, 011506 (2005). 

[40] S. S. Dukhin, "Electrokinetic phenomena of the second kind and their applica- 
tions". Adv. Coll. Interf. Sci. 35, 173 (1991). 

[41] S. S. Dukhin and N. A. Mishchuk, "Disappearance of limiting current phe- 
nomenon in the case of a granule of an ion-exchanger", [in Russian] KoUoidn. 
Zh., 51, 570 (1990). 

[42] S. S. Dukhin, N. A. Mishchuk, and P. B. Takhistov, "Electroosmosis of the sec- 
ond kind and unrestricted current increase in the mixed monolayer of an ion- 
exchanger", [in Russian] Kolloidn. Zh., 51, 540 (1989). 

[43] S. S. Dukhin and B. V. Derjaguin. Electrophoresis, 2nd edn. [in Russian] Nauka, 
Moscow, 1976. 

[44] E. K. Zholkovskii, M. A. Vorotynsev and E. Staude, "Electrokinetic instability 
of solution in a plane-parallel electrochemical cell", J. Membrane. Sci. 181, 28 
(1996). 

[45] M. Z. Bazant and T. M. Squires, "Induced-charge electro-kinetic phenomena: 
Theory and microfluidic applications", Phys. Rev. Lett. 92, 066101 (2004). 

[46] I. Rubinstein and B. Zaltzman, "Electro-osmotically induced convection at a 
permselective membrane, " Phys. Rev. E 62, 2238 (2000). 

[47] I. Rubinstein and B. Zaltzman, "Electro-osmotic slip of the second kind and 
instability in concentration polarization at electrodialysis membranes, " Math. 
Mod. Meth.Appl. Sci., 11, 263 (2001). 

[48] I. Rubinstein and B. Zaltzman, "Wave number selection in a nonequilibrium 
electroosmotic instability", Phys. Rev. E 68, 032501 (2003). 



40 



[49] I. Rubinstein, B. Zaltzman, and I. Lerman, "Electroconvective instability in con- 
centration polarization and nonequilibrium electro-osmotic slip", Phys. Rev. 72, 
011505 (2005). 

[50] T. Pundik, 1. Rubinstein and B. Zaltzman, "Bulk electroconvection in elec- 
trolyte", Phys. Rev. E., 72, 061502 (2005). 

[51] V. I. Zabolotsky and V. V. Nikonenko. Ion transport in membranes, [in Russian] 
Nauka, Moscow, 1996. 

[52] L. I. Sedov. Similarity and dimensional metods in mechanics, 4-th edn. Academic 
Press, New York, 1977. 

[53] Yu. M. Shtemler, "Stability of unsteady viscous flows", [in Russian] Izv. A. N. 
USSR. Mech. Zhidk. i Gaza. 4, 138 (1981). 

[54] C. C. Lin. The theory of hydrodynamic instability, Cambridge University Press, 
Cambrige, 1967. 



41 



Figure Captions 

Figure 1 Schematic V-C curve for ion-exchange membranes. A, B and C represent under- 
hmiting, limiting and over-limiting regimes, respectively; * depicts the threshold 
of the over-limiting regimes; Cj and Cn are the regions of regular and irregular 
current oscillations; j and jum are the current density at the bottom membrane and 
its limiting value; AV^ is the potential drop. 

Figure 2 Schematic for limiting processes in a steady-state equilibrium: a) electric field 
E = d^/dy; b) volume charge density p= — c^; c) ion concentrations c = c'^ and 
c = c^. The dashed line is the outer solution in the small vicinities of y = (the 
double ion boundary layer) and y = ym (the boundary layer between the volume 
space charge region and the diffusion layer). The solid and dashed lines depict 
exact and asymptotic solutions for small e. 

Figure 3 Marginal stability curves at different e, numerics. The under-limiting regimes are 
stable (adopted from Demekhin et al. jTU]). 

Figure 4 Self-similar charge p = c^—c~ (a) and concentration c"*" distribution (b) for different 
e. Curves correspond to the following values of electric potential drop: 1. A$ = 50; 
2. A<l> = 100; 3. A<l> = 150; 4. A$ = 200. 

Figure 5 VC curves. Comparison of the self-similar asymptotics (dashed lines) and the exact 
numerics (solid lines) for several e. 

Figure 6 Distribution of the charge density p in space for several moments of time, AV = 100 
and z/ = A/)/I/ = 0.001. Solid line corresponds to the numerical solution of the 
unsteady problem (I32i) - (l36|) . Two-parametric self-similar solution of (128|) - (13T]) 
with e and A<l> parametrically changing with time is depicted by dashed lines. 

Figure 7 Typical evolution of 1 / j{t). I is the short-time region of influence of the initial data; 

II is the intermediately-long-time region of self-similarity, and III is the region of 
influence of the upper membrane 1/=1, (AK = 100 and z/ = A^j/L = 0.001). 
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Figure 8 VC curves. Solid line relates to flMj) . while dashed line - to the simplified version 



Figure 9 Negative ion concentration (a), c , and electric field (b), if, vs y for several mo- 
ments of time. Comparison of numerics f|T3|) - f|T8l) (solid lines) with self-similar 
asymptotics (|19D-([50D (dashed lines) for (a) Ay = 50, z/ = 0.001: 1. t = 0.0001; 2. 
t = 0.001; 3. t = 0.004; 4. t = 0.015; 5. it = 0.5; and (b) AV=100, i/ = 0.0005: 1. 
t = 0.001; 2. t = 0.003; 3. t = 0.03. 

Figure 10 Comparison of the V-C characteristics for self-similar solution and numerics, = 
50, V = 0.00005. Numerics flT3|) - f|T8|) is shown by squares and the universal self- 
similar VC curve - by solid line. The region of the solution self-similarity 
ti< t <t2, ti 3> z^^/4, t2 <^ ts = L^/(4Z)). Arrows depict the time growth along 
the voltage-current curve. 

Figure 11 Numerical points for several values of the potential drop between the membranes, 
AV^, shrink into the universal VC curve J{AF) given by flMl) : a) z/ = 0.001, b) 
i/ = 0.0005. 

Figure 12 Components of scaled slip velocities u and v. 

Figure 13 Stream-lines for A$ = 100, e = 0.001 and different values of a. 

Figure 14 Marginal stability curves for self-similar asymptotic in e solution fll06p - fll07p 
(dashed lines) and exact numerical solution (solid lines) [TO] . 

Figure 15 Marginal stability curves in 1/e and a coordinates. 

Figure 16 The universal functions a) G = G{r), b) a = a{r). 
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Table Caption 

Table Critical parameters which separate stable and unstable regions (tangent ko, wave 
number and l/^o); and Q and R are functions of A$. 
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A. Electrostatic solution in the space charge region 

Substituting a new independent variable ^ = a {rjm — v) i^^to (PT]) - (|^ results in: 

F = for ^ = (114) 

F = for i = r. (115) 

where r = ar/^. The additional conditions which originate from (192 p are: 

2a 1 dF 

T^'de^ at e^O, (116) 



, J d^F dF , , , , 

aJ- — 2r-— + — =-J for ^ = r. (117) 



2Vm \ de di 
We present the solution of f lll3p - flll5l) as 

F = B^fiO (118) 
where B is constant. The equation ( I113P turns into 

er + If" - e/' - ^/ = 0. 

The general regular solution of this equation, up to a constant factor, has the form 
F = -By^ (/i (^) + a/2 (0); where /i and /2 are two linearly independent solutions of the 
equation, a is determined below. Solutions /i, /2, can be easily found numerically. In- 
stead, for comparison convenience with long-wave asymptotics, the solution is presented 
as power expansions 

00 00 

/i(0 = E«2ne^ /2(o = E«2n+ie'"+\ (119) 

n=0 n=0 

where ~ ^/k\, and the series have infinite radius of convergence 



_ 1 A4A;-1 



1 -rT 4fc + l _ 

Note that the long wave limit is described by a finite term expansion (11190 . 
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(120) 





Figure 16: The universal functions a) G = G{r), b) a = a{r) 



Assuming 



/ = /i + a/2 = 



at some arjm = r (a), the second BC flllSp for F is satisfied. Let us invert the function 
r = r(a), 

f. (r) ^ ^ 

121 



a[r) = — J. f ^ , r = arjn 



f2{r) 



Thus we obtain, taking into account flllSp . that the function F satisfies the first two 
boundary conditions ( I114p - flll5p . The constant B is obtained by substitution of f lllSp 
into ffTTTD . 



Vn 



G(r) J 



G{r) 



1 



(122) 



a J' ^' ' 2rf"{r) + 3f'{r)' 
where /', /" are taken at ^ = r and are functions of r, G = G (r) is the universal function 
shown along with a(r) in Fig. [T6l 

Then we get the unknown function F, 



2J 
V a 



(123) 



Substituting f ll23p into (!93|) . we obtain the relation for the outer expansion of the charge 
density /), 



p = eV2Ja^ fjr, 
which has a singularity at r] = rjm- 



(124) 
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Furthermore, from flll9p -f lT2T|) . the asymptotic expansions for a(r) and G{r) as r — )■ 
can be found 

«W = - Sr^''^' 2n^V (125) 

For small r, which corresponds to the long-wave approximation {a — )■ 0),we obtain, 

/ N 1 19 169 o 1 ^ 

a(r) = rH r"* H ~ — , r 0. (126) 

^ ^ r 105 24255 r' ^ ^ 

From (11221) . in the first approximation with respect to r, 

1 1 59 1 

G(r) = -r H H ~ - r, r 0. (127) 

^ ^ 3 210 415800 3 ' ^ ' 



and 



Approximation ffTTTjl -f ing]) is valid at r < 1.5 ^ 2. 

The asymptotics of /i and /2 at ,^ — )■ oo are useful in the short-wave approximation 
a — 7- oo of the stability problem. They can be readily found by asymptotic expansions of 
the equation (I113p : 

„ consti . . const2 f . 
/i ~ ^ eS /2 7^ ' 



At r — 7- oo 



where 



a(r)~a„ G(r)~c,v^, (129) 



a* ^ -0.72948, c. ^ 0.66468. 



B. Components of the slip velocity 

Now the obtained function F can be used to solve the hydrodynamic problem (l94l) . 
flUS]) . f lUB]) and find the stream function perturbation ^ . The slip velocity components 
^ ('7m) = "^Kre and ^ (//„) = f/^ specify the boundary conditions for the problem in the 
electro-neutral region rjm < t] < oo. 

Taking the variable ^ = ^(r^rn — v) with r = arj^, and utilizing relations (I122p . (11231) . 
the problem fl94|) - (l96|) can be reformulated as 
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where 



^(e,r) = ^^, C = 4^G(r)J. (132) 



Vm and Um are used as the boundary condition for the hydrodynamic problem in the 
electro-neutral region t] > rj^, 

= Cv (r) = Ci) (0, r) , f/„ = -%aCu (r) = -^aC^ (0, r) . (133) 

The function / in the right hand side of fll30p is 

, fd'f A ^ df def 



where 



oo 

„ _ 7, f 2n+l , _ ^2n 

(2n + 1) (4n + 5) 



n=0 

oo 



yi ^ 2n^ , in 2n(4n + 3)' 



(134) 



n=0 

Then the solution of (11301) — ( ITM]) can be presented in the form 

^ = (A + B^e^ + , = ^1 + a^2 , (135) 

where and ^e^ — are two linearly independent solutions which satisfy the first two 
conditions (I13ip at ^ = 0, A, B — some constants, ipp is a particular solution of f ll30p 
which satisfies zero conditions 

Functions ipi and V'2 ~ are the odd and even components of ipp, which can be found 
from solutiion (I130p with right-hand side gi and g2 correspondingly and boundary con- 
ditions (I136p . They can be presented as expansions 

n-2 / 

[n — , 

(2n + 1)! 

(137) 

V- .2n ^in-k-l){2k)\ ^ 

n=2 k=0 ^ ^' 
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Wl- C2n+lt, , C2„+l - 02fc+l, 



n=2 fc=0 
oo n—2 



The constants A, B are found from the BC's at ^ = r, 

A = - (1 + r)i;,{r)\ e-^ B = (^,(r) - ^ ) e'^ (138) 



Taking into account (11331) yields 



" " V " ^ ""^'^"^'^V ^^^^^ 

u = A + B=(^{r-l)^^^i^-rijp{r)^e-', (140) 

The obtained relations provide a complete description for the slip velocity components 
Um and Vm according to (I133p . 

Finally, let us present the Taylor series near r = for the slip velocity in the long-wave 
approximation. Using (11261) . (11351) and (11371) . we can write the following relation, 

where the series coefficients are calculated according to (11201) and (I134p . (I135p . Expanding 
the exponents in (I139p . (I140p in series in the vicinity of r = 0, re-expanding them 
and differentiating the above relation with respect to r, the functions u{r) and v{r) are 
obtained, 

1 3 1 4 27 5 

vir) = — r H — r r + ■ ■ ■ , 

8 6 200 .142) 

uir) = —r r H r + • • • . 

^ ^ 6 3 504 

All the series in (I14ip converge at any r, while the series (I142p have finite convergence 
radius because of above mentioned hidden singularity of the universal function a = a(r) 
in the complex plane. 

We find the limiting values of u and v, 

v-^Voo^ -0.27245, M ^ Moo ~ -0.12437 as r ^ oo. (143) 

Moreover, u changes sign at r = arjm ~ 1.59523. The dependence of u and f on r is 
shown in Fig. [121 At r — )■ results can be obtained in a closed form; they coincide with 
obtained in |i7] . 
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